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Abstract: Neutralinos arise as natural dark matter candidates in many supersymmetric 
extensions of the Standard Model. We present a novel calculation of the neutralino relic 
abundance in which we include all so called coannihilation processes between neutralinos, 
charginos and sfermions, and, at the same time, we apply the state of the art technique to 
trace the freeze-out of a species in the early Universe. As a first application, we discuss 
here results valid in the mSUGRA framework; we enlight general trends as well as perform 
a detailed study of the neutralino relic densities in the mSUGRA parameter space. The 
emerging picture is fair agreement with previous analyses in the same framework, however 
we have the power to discuss it in many more details than previously done. E.g., we find 
that the cosmological bound on the neutralino mass is pushed up to ~ 565 GeV in the stau 
coannihilation region and to ~ 1500 GeV in the chargino coannihilation region. 



Keywords: 



supersymmetry, dark matter. 



Contents 

[l]. Introduction |] 

^ Supersymmetric model ||| 

^ Relic density calculation |B] 

|3.1| The density evolution equation and thermal averaging H 

3.2 A few examples of coannihilation effects ^ 



The role of coannihilations in the mSUGRA framework 13 

|4.1| Slepton coannihilations 13 

[4.2| Chargino coannihilations 19 

Stop coannihilations 21 



|5). Conclusions 24 

|6]. Acknowledgements 24 

|A] . Included coannihilations 26 

|B|. A note about internal degrees of freedom 28 

Example models 29 



1. Introduction 

The latest years will be remembered in the history of Science as those that marked the 
entrance into the era of precision Cosmology. A number of experiments have been pin- 
ning down the values of cosmological parameters to a level of precision hardly foreseeable 
just a decade ago, with perspectives from upcoming measurements even more spectacular. 
Most remarkably, experiments with different focus, as well as exploiting complementary 
techniques, have all collected data which point to one single overall-consistent picture, a 
"concordance" model ]|] in which the Universe is flat with about 30% of its present av- 
erage energy density in matter and about 70% in some form with negative pressure (a 
cosmological constant or dark energy). More precisely, from the combined analysis of the 
latest data on the cosmic microwave background and large scale galaxy surveys, the cold 
dark matter (CDM) and baryonic contributions have been recently estimated [|| to be 
V-CDMh 2 = 0.115 ± 0.009 and Q b h 2 = 0.022 ± 0.002, respectively (here Q is the ratio be- 
tween mean density and critical density p c = 1.879 x 10 _29 /i 2 g/cm 3 , and h is the Hubble 
constant in units of 100 km s _1 Mpc -1 ; the analysis in Q gives h = 0.665 ± 0.047). 
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Among the ideas which have been put forward to account for the CDM term, the most 
natural solution is probably the scheme in which CDM appears as a thermal leftover from 
the early Universe: in this context, stable weakly interacting massive particles (WIMPs) 
are ideal CDM candidates, as their thermal relic abundance is naturally of the order of the 
measured one. Given the accuracy of the current and future measurements of Qcdm, it 
is useful to have an equally accurate calculation of the relic abundance of WIMPs. This 
paper continues our program of accurate relic density computations. 

Since the late seventies, when the idea of WIMP dark matter was formulated |, |,§, 
the computation of the WIMP relic abundance has been constantly refined. One important 
step was to recognize the role of coannihilation effects || 0] : if in the particle physics theory 
one considers, a stable WIMP appears together with a slightly heavier particle into which 
it can transform, when computing the present density of the lightest particle one needs to 
retrace the thermal history of both particles simultaneously. Such an effect is common even 
for the most popular WIMP dark matter candidate, the lightest supersymmetric particle 
(LSP) in supersymmetric extensions of the particle physics Standard Model. When the 
LSP is the lightest neutralino, the mass splitting between the LSP and the next-to-lightest 
supersymmetric particle may in some cases be small. Coannihilations are then important. 
Their inclusion in the computation of the neutralino relic density has been the subject 
of numerous studies, which vary in the degrees of refinement of the method implemented 
to compute the relic abundance and in the selection of which particles to include in the 
coannihilating set. 

The novel calculation that we present in this paper is included in a new extended 
version of the DarkSUSY package [||, which will be publicly available in the near future. 
DarkSUSY now allows for the most generic coannihilation effect in the framework of the 
minimal supersymmetric extension of the Standard Model (MSSM), at the same time 
applying the state of the art technique to trace the freeze-out of a species in the early 
Universe Q. In fact, on one side, we fully solve the density evolution equation (which 
determines the evolution of the number density of neutralinos) numerically, including all 
possible resonance and threshold effects, and avoiding approximations in the thermally 
averaged cross sections (such as the expansions in terms of the relative velocity that is often 
applied). On the other side, extending the work of two of us [1C] who first applied such 
a technique to the case of coannihilations with charginos and neutralinos, we include here 
the possibility of having coannihilations with all sfermions as well. As a result, assuming 
masses, widths and couplings of particles in the MSSM are given with an adequate precision, 
we provide here a tool to compute neutralino relic abundances with an estimated precision 
of 1% or better. 

Compared to other recent calculations, we believe this is the most accurate calculation 
available at present. The standard lore so far has been to calculate the thermal average 
of the annihilation cross section by expanding to first power in temperature over mass 
and implementing an approximate solution to the evolution equation which estimates the 
freeze out temperature without fully solving the equation (see, e.g., Kolb and Turner 
Sometimes this is refined by including resonances and threshold corrections M. Among 



recent studies, this approach is taken in e.g. Refs. [Jl^, 12]. Other refinements include 
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e.g., solving the density evolution equation numerically but still using an approximation to 
thermal effects in the cross section [ 14 , [15], [If], [l?], |l8| , or calculating the thermal average 
accurately but using an approximate solution to the density equation [|^, At 
the same time, only in a few of the quoted papers the full set of initial states has been 
included. As already mentioned, the present calculation includes all initial states, performs 
an accurate thermal average and gives a very accurate solution to the evolution equation. 
Though the inclusion of initial state sfermions in the DarkSUSY package is a new feature 
introduced in the present work, other groups [22, 23, 24 1 have earlier introduced some 
sfermion coannihilations in an interface with the old DarkSUSY version. 

The DarkSUSY package has been written in a very general and flexible format, so that 
it can be used for any theory embedded in the MSSM. As a first example, we will present 
here results valid in the minimal supergravity (mSUGRA) framework. The mSUGRA 
framework is the framework considered in most of the previous analyses. Although rather 
restrictive in the way parameters are set, it is sufficient to enlight most coannihilation 
effects which can emerge in the MSSM. 

The outline of this paper is the following: in the next section we discuss the supersym- 
metric model we work in; in section |3| we review the framework in which we calculate the 
relic density including coannihilations. In section 4 we examine the effects of coannihilation 
in detail, stressing the physical insights of the results. 



2. Supersymmetric model 

The particle physics model implemented in the DarkSUSY package is a minimal supersym- 
metric extension of the Standard Model, built with N = 1 generator of supersymmetry 
and containing the smallest possible number of fields. We restrict ourselves to the case 
in which the LSP is the lightest neutralino, i.e. the lightest of the four mass eigenstates 
obtained from the superposition of the supersymmetric partners of the neutral gauge and 
Higgs bosons. The supersymmetric part of the spectrum contains also two chargino mass 
eigenstates, the gluino and the scalar superpartners of leptons and quarks. The Higgs sec- 
tor needs two Higgs doublets, which, after electroweak symmetry breaking, give five scalar 
fields, denoted as H^,!!®, H% an d H® (for a more detailed description of the model and of 
our conventions see f^ ). 1 

As needed in the relic density calculation, we perform the computation of all two-body 
final state cross sections at tree level for all initial states involving neutralinos, charginos, 
sleptons and squarks. We do not neglect mass terms for fermions in final states (as some- 
times done in the literature), nor implement any expansion to first order in relative mo- 
mentum of the incoming particles (as most often done in the past). 2 A list of all processes 
included is given in Appendix |A|. 

1 Apart from sfermion coannihilations, another major improvement in the DarkSUSY version that we 
use for this work is the inclusion of one-loop corrections to the Higgs widths, with formulas taken from 
Refs. H 0. 

2 In the current version of DarkSUSY we do not include processes with the gluino in the initial state or 
any flavour changing processes in the sfermion coannihilations; these might be included in future work but 
we do not expect any result in the present analysis to be affected. 
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We compute analytic expressions for the amplitudes using standard Feynman rules with 
generic expressions for the vertex couplings (e.g. ig^^Pi + igjpjjr,Pji for a vertex involving 
two fermions and a scalar) . The amplitudes for all the different types of Feynman diagrams 
are obtained with the help of symbolic manipulation programs. For initial states involving 
neutralinos and charginos, helicity amplitudes for each type of diagram are computed 
analytically with Reduce p8| (this is described in [10] 3 ). For the diagrams including 



sfermions in the initial state, Form [£9| is used to analytically calculate the amplitudes 
squared. The output from Form is then converted into Fortran with a Perl script. All 
of these analytic formulas can be converted into a compact form, but we do not consider 
it useful to reproduce those expressions here. 

The actual values of the MSSM vertex couplings are introduced during the numerical 
calculation. We recomputed all vertex couplings from the MSSM lagrangian, and checked 
them against standard literature (e.g. |3(| |3l]]; DarkSUSY also contains couplings which do 
not appear in the literature, namely those involving generic intergenerational mixing in the 
sfermion sector). 

When giving explicit examples of computations of the neutralino relic abundance 
we restrict ourselves in this paper to mSUGRA models. Under the assumption of uni- 
versality at the grand unification scale, the mSUGRA action has five free parameters, 
mi/2,Tno,sign(fJL),AQ and tan/3. The parameters muz, mo and Aq are the GUT unifi- 
cation values of the soft supersymmetry breaking fermionic mass parameters, scalar mass 
parameters and trilinear scalar coupling parameters, respectively (of the trilinear couplings, 
only At, A}, and A T differ from zero). The absolute value of the Higgs superfield param- 
eter n follows from electroweak symmetry breaking, but its sign is free. Finally, in the 
Higgs sector, tan/3 denotes the ratio, v^/vx, of the vacuum expectation values of the two 
neutral components of the SU(2) Higgs doublets. Our convention on the sign of \i is that 
// appears with a minus sign in the superpotential (i.e. following the convention used in 
e.g. p0[| ), while the definition of sign and dimension for the Aq's are such that At appears 
in the stop mass squared matrix as off-diagonal elements of the form (At — //cot /3)mt (and 
analogously for A^ and A T ). 

GUT scale values of the soft breaking parameters, as well as of the gauge and Yukawa 
couplings, have to be evolved down to the weak scale. To do that we make use of the 



ISASUGRA RGE package in the ISAJET 7.64 software || 4 . The interface to DarkSUSY 
is such that the whole SUSY spectrum of masses and mixings as given in the output at the 
weak scale in ISASUGRA is given as input in the DarkSUSY MSSM spectrum. Regarding 
Standard Model masses and couplings, we have implemented 1-loop renormalization group 
equations. We fix the top pole mass at 174.3 GeV, which is the central value stated by the 
Particle Data Group 2002 |33|. For the c and b quarks we input the running quark masses 

3 Compared to the analysis in (To[ , we have here corrected a sign error in the chargino - (down-type) 
fermion - (up-type) sfermion vertex with clashing arrows. Consequently, the discrepancy between Dark- 
SUSY and micrOMEGAs (model B in has now disappeared. 

4 The ISASUGRA output is rather unstable in regions where the convergence is slow (e.g. in the focus 
point region where the lightest chargino is nearly degenerate with the lightest neutralino); to improve 
the stability, we have converted ISASUGRA to double precision and increased the requirements on the 
convergence. 
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m q (m q ) and choose the values of 1.26 GeV and 4.2 GeV respectively [33, 34]. To a good 
approximation, the RGE scale at which (co) annihilations occur is twice the LSP mass, and 
therefore we evaluate the gauge and Yukawa couplings at this scale. 

We have chosen to use the ISASUGRA RGE code since it is a widely used program 
that is kind of a standard in the field. However, there are other RGE codes on the market 
(e.g. SOFTSUSY ||, SPHENO §| and SUSPECT @), all of which solve the RGEs 
with a different level of sophistication. We do not want to go into the details of these 
different packages here, but instead refer the interested reader to a recent comparison in 
Ref. |3q] . The main message we want to convey here is that different RGE codes give 
slightly different results. E.g. the sparticle masses typically differ by a few GeV between 
the different codes. 

When scanning over the parameter space, we compare with accelerator constraints, 
implementing limits recommended by the Particle Data Group 2002 (PDG) p3| . We have 
not included limits which are still preliminary as of this writing (most important among 
the preliminary bounds is the limit of about 104 GeV on the chargino mass; we have 
nevertheless indicated its effect in some of the figures). For the Higgs's H® and H® we 
have implemented the tan/3 dependent mass limits as provided by the PDG in the most 
conservative setup. This is particularly important for the mass of since for a substantial 
tan/3 interval this constraint is less restrictive than the lower mass limit (114.3 GeV) on 
the Standard Model Higgs. In this paper, which focusses on coannihilations, we do not 
emphasize the role of the b — > 57 constraint (we refer the interested reader to recent 
dedicated analyses in the same mSUGRA context, e.g., Refs. p9|). 



3. Relic density calculation 

The importance of coannihilations in the computation of the density of a relic particle was 
recognized by Binetruy, Girardi and Salati ||], and independently by Griest and Seckel Q. 
In Edsjo and Gondolo [10], this was further analysed and put in a form that allows for an 
accurate treatment of coannihilations in the same basic framework as without coannihila- 
tions. We do not repeat every step of that calculation here. Instead, we only briefly review 
the results we need for this paper. For more details, we refer the reader to Ref. flcfl . 

3.1 The density evolution equation and thermal averaging 

We consider the coannihilation of N species of particles (xi, i = 1, ■ ■ ■ , N) with masses rrii 
and internal degrees of freedom (statistical weights) gi (see Appendix [B] for conventions 
on degrees of freedom adopted in this paper). We order the masses in increasing order, 
m\ < ni2 < • • • < niN, and use mi and m x interchangeably for the mass of the lightest 
particle. If the lightest particle is stable and the others decay into it, instead of considering 
the thermal history of each particle separately, we can follow the evolution of the sum of 
the number densities. The problem is then formulated in terms of the density evolution 
equation [Q, |o| 

— = -3Hn - (a cS v) (n 2 - n^) (3.1) 
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where ((T^v) is the effective thermally-averaged annihilation cross section, H is the Hubble 
parameter, and n is the total number density summed over all coannihilating particles. The 
effective thermally-averaged annihilation cross section is 

WefiV) = 4~ > ( 3 ' 2 ) 

where n eq is the equilibrium number density, which in the Maxwell-Boltzmann approxima- 
tion (which is a good approximation for our case) reads 

T \ -« 2 f m i\ 

neq = 2^2 2^ 9iTn i K2 \~r ) ' ^ 3 ' 3 ^ 

i 

and A is the annihilation rate per unit volume at temperature T 

a _ 9i T f°° , 2 ttt _ ts ( Vs s 



4vr 4 



J^dpeSplgWesKj^Pj . (3.4) 



Here Ki(x), (i = 1,2), are the modified Bessel functions of the second kind of order i; 
s = 4pg ff + 4m 2 is the usual Mandelstam variable giving the center-of-mass energy squared 
for the XiXi system; p e g is the center-of-mass momentum for the XiXi system; and W e s is 
the effective annihilation rate obtained by summing over all annihilation and coannihilation 
channels, 



Vij 9i9j w \- / [s - (mi - mj) 2 ] [s - (m, + m 3 -) 2 ] cfegj 



^ Pll 9i ^ V s(s - 4mf ) 



For the coannihilation of particles i and j, Wij is the annihilation rate per unit volume and 
unit time given by 5 

Wij = 4pij\fs(iij = AaijJ (pi ■ pj) 2 — mfm| = AEiEjaijVij. (3.6) 

where 

[s - (m,j + rnj) 2 ] 1/2 [a - (rm - rrij) 2 ] 1/2 

is the common magnitude of the 3-momentum of particle i and j in the center-of-mass 
frame of the pair XiXj- Defining Wij(s) = for s < (rrtj + rrij) 2 , the radicand in Eq. (|3.Ej ) 
is never negative. The effective thermally averaged annihilation cross section can then be 
written as 

(vezv) = (3.8) 

^ gim* 2 \ T ) 



This expression is very similar to the one in the case of no coannihilations given in Gondolo 
and Gelmini || (and correctly reduces to it in the absence of coannihilations). The only 



5 The quantity Wij in Ref. |40| is Wij /A, and is therefore one-fourth of the annihilation rate per unit 
volume and unit time. 
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differences are in the denominator and in the replacement of the annihilation rate with 
the effective annihilation rate. The key feature in Eq. ( |3.8| ) is the definition of an effective 
annihilation rate independent of temperature, with p e R as integration variable. This gives 
a remarkable calculational advantage, as W e fj can be tabulated in advance, before taking 
the thermal average and solving the density evolution equation. 

The steps we implement to compute the relic density are the following. We first 
rephrase the density evolution equation ( |3.1[ ) as an equation for Y = n/s, with s being the 
entropy density. Then we tabulate W e s including thresholds, resonances and coannihila- 
tions, and spline it. As a third step we solve the density evolution equation starting from 
the boundary condition which sets particles in equilibrium at the temperature T = rn x /2 
(a specially-devised implicit method is used in the numerical integration of the density 
equation). This implies actually a double integration since at each temperature step we 
need to calculate the thermal average (cr e ffv) (this integration is relatively fast as we use 
the W e g tabulation). The integration of the density equation is performed until Y has 
reached a constant value (the freeze out value) which we finally convert into the value of 
the relic abundance. 

As stated earlier, we estimate our calculation of the relic density to be accurate to 
about 1%. We base this estimate on the precision with which we do the tabulation of the 
effective annihilation rate W e g and on the precision with which we perform the numerical 
integrations. For example, we explicitly make sure that resonances and thresholds are 
tabulated and integrated with such a precision that the end result is accurate to at least 
1%. One should keep in mind that this is the accuracy of the calculation of the neutralino 
relic density starting from given values of masses, widths and couplings of MSSM particles. 
In the mSUGRA framework considered here, the accuracy of the sparticle masses for a given 



set of input parameters is less than 1% (see e.g. |38(] for a comparison of different RGE- 
codes) and the main uncertainty in evaluating the neutralino relic density in mSUGRA 
thus comes from the RGE code ISASUGRA. 

3.2 A few examples of coannihilation effects 

For illustrative purposes only, we rewrite Eq. ( |3.8| ) in the form: 

(a cS v ) = f dp cS W f^ cf£ ^ K(p cS ,T) , (3.9) 
Jo 4£ eff 



where E e g = yP^g + m x * s ^ e ener gy P er particle in the center of mass for xiXi anni- 
hilations. The term we factor out, W c g/AE^ S , can be thought of as an effective av term 
(compare with Eq. (|3.6|)). In thep e ff — * limit, W e g/4E^g reduces to the XiXi annihilation 
rate at zero temperature, which is the relevant quantity in indirect dark matter detection 
searches. The term k we introduced in Eq. ( |3.G| ) contains the Boltzmann factor and the 
phase-space integrand term and can be regarded as a weight function, at the temperature 
T, that selects which range of p e g is important in the thermal average. As the phase-space 
integrand term dominates at small p e g and makes k go to in the p e g — * limit, k shows a 
peak at an intermediate p e g and then rapidly decreases due to the Boltzmann suppression; 
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Figure 1: The effective annihilation cross section a) with coannihilations and b) without coan- 
nihilations for model A (specified in Table |2| in Appendix |C]). The solid line shows the effective 
annihilation cross section W c s/AE^ S as a function of momentum p c g , while the dashed line shows 
the thermal weight factor K(p e g,T). The thermally-averaged annihilation cross section is the inte- 
gral over p e g of the product of the two. Note that when including coannihilations, not only new 
thresholds appear, but the freeze-out temperature is also changing, meaning that we sample a dif- 
ferent region of the annihilation cross section. For this model, the relic density with coannihilations 
is fi x .coann^ 2 = 0.135 and that without is fl x . no coann /! 2 = 1.43. 

the position and height of the peak depends on the temperature considered and on the 
particles involved. 

We are now ready to show some examples of coannihilation effects. As already men- 
tioned, the examples we display have the lightest neutralino as the LSP and are in the 
mSUGRA framework. In Fig. |l]a we consider a case in which the neutralino, with mass of 
about 400 GeV, is nearly mass degenerate with the lightest stau. The lightest selectron, the 
lightest smuon and the lightest stop are relatively close in mass as well. (To fully specify 
the example models we present, the model parameters and some properties are given in 
Table ^ in Appendix [C]. The model in Fig. [I] is model A in that table.) The solid curve 
shows W e ff/4Eg ff , and one can nicely see coannihilations appearing as thresholds at 
equal to the sum of the masses of the coannihilating particles (just as final state thresh- 
olds do). As usually happens when considering coannihilation effects with neutralinos as 
the LSP, the x?~Xi contribution to W e R is small compared with the one provided by the 
coannihilating particles. The role of coannihilating particles can be quantified better with 
a look at the function k (dashed curve, in units of GeV -1 , and with relative scale shown 
on the right-hand side of the figure). The factor n is plotted at the freeze out temperature, 
defined as the temperature at which the abundance of the relic species is 50% higher than 
the equilibrium value 6 , in this case T = m x /24.3. On the top of the panel, the tick mark 

6 This is given here for illustrative purposes only; it is never actually exploited in the full computation 
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labelled '0' indicates the position of the momentum p^ x corresponding to the maximum of 

(n) — 

k, while the other tick marks indicate the momenta p eS at which k is 10 n of its maximum 
value, k(p^)/k(p^ x ) = 10~ n . The tick marks provide a visual guide to the interval in p c g 
which is relevant in the thermal averaging. The integral of the product of W e d / 4:E 2 S and k 
gives (cr e fju ) thermally averaged at the freeze out temperature (shown in the figure as a hor- 
izontal dotted line) . This is the quantity which is sufficient to get a rough indication of the 
neutralino relic abundance through the rule of thumb [^|] Q, x h 2 ~ 1CP 27 cm 3 s^ 1 / ((T c gv) . 

In Fig. |l]b we consider the same model but ignore coannihilation effects. One can 
see that W e s/4:E 2 S is now, on average, much smaller, and therefore one can expect the 
relic abundance to be higher. This is indeed the case, with a shift from £l x h 2 = 0.135 
including coannihilations (left panel) to Q x h 2 = 1.43 when coannihilations are neglected 
(right panel). Note, however, that the change in Q. x h 2 is smaller than what one would 
naively expect from comparing the solid curves in the two panels. This is due to the fact 
that there is a significant change in the freeze out temperature as well, from T = m x /24.3 
to T = m x /21.7. The weight function k for this new temperature is shown in the figure, 
and comparing it to the one in the left panel, one clearly sees the change in normalization 
(partially due to the change in the number of degrees of freedom involved in the two cases, 
see the denominator in Eq. Q3.8| )) and in width (for the latter effect, note the shift in the 
scale shown on the top of the figure, while the displayed range in p e s has been kept fixed). 
The net result is that (a c sv) at the freeze out temperature is lowered by just about an order 
of magnitude, and then the increase in the relic abundance is of the same order. From this 
discussion it is evident that a very accurate solution of the density evolution equation is 
needed to claim good accuracy on the estimate of the relic density. 

We mentioned that usually, when considering neutralino dark matter, W e s increases 
sharply when coannihilating particles are included. The reason is that the coannihilating 
particles typically have non zero electric or colour charges, while the neutralino interacts 
only weakly. In the MSSM there are a few exceptions to this general trend of increasing W e s 



(see [10]), and we find one such exception in the mSUGRA framework as well. In Fig. ga we 
show a model (model B in Table || in Appendix |C]) for which the neutralino mass is slightly 
below half of the H® mass and tan j3 is large. The Xi~Xi annihilation rate is dominated by 
the s-channel H® resonance and is quite large. (There is also a H® resonance at this mass, 
but it is subdominant with respect to the #3 resonance). The effect of stau coannihilating 
particles comes on top of that, but the total contribution to W e s is just about of the same 
order as from the x?"X? term. If we now go to the case when coannihilations are neglected, 
Fig. Qb, we see that the freeze out temperature remains about the same, but still there is 
a shift in the normalization of the weight function due to different numbers of degrees of 
freedom in the two cases (see the denominator in Eq. (|3.8|)). We find that the net effect 
is a slight increase in (cr e ffv), which suggests that the relic abundance should be smaller if 
coannihilations are neglected. This is confirmed by the full solution of the density evolution 
equation: the calculation that includes coannihilations yields Q x h 2 = 0.155, while if stau 
coannihilation is neglected one obtains Q x h 2 = 0.137. Therefore one should keep in mind 



since we solve the density evolution equation numerically. 



- 9 - 



10 



-23 -1 



-1 -2 -3 -4 -5 -6 -7 -8 



F i i i i i i i r 

With coannihilations 
H? res. 



10 



-28 - 
10 r 



10 




T = m/24.3 
nil =0.155 
m* = 646.5 GeV 
m== 665.1 GeV 
m„, = 1323.1 GeV 



Id' JT 10 



-23 -1 



-1 -2 -3 -4 -5 -6 -7 -8 



100 200 300 400 500 600 700 
P eff [GeV] 



F i i r 

Without coannihilations 
H? res. 



10- v. l£ 10 



- io- 4 



10 



-28 - 
10 r 



10 




T = m /24.0 
£2h =0.137 
m: = 646.5 GeV 



10 1 j-i 

U 

10" 2 iT 



- 10'" 



100 200 300 400 500 600 700 
P e „ [GeV] 



Figure 2: We here show the effective annihilation cross section versus p e ff for an example (model 
B in Table ^ in Appendix ^]), where coannihilations increase the relic density (in this case from 
0.137 to 0.155). 

that although the general trend is that co-annihilation effects lower the neutralino relic 
abundances, in some cases they can increase it. 

Note that due to the Boltzmann suppression of heavier particles, we do not need to 
include all supersymmetric particles in the calculation. By extensive scans of the mSUGRA 
parameter space and estimating the effects of including different particles in the calculation, 
we have found that a convenient criterion to have an accuracy of 1% or better on the relic 
density in cases that are cosmologically interesting is to include all supersymmetric particles 
with a mass below 1.5m x . This is shown in Fig. || where we plot the relative difference in 
relic density without coannihilations and with coannihilations, 

^Xi nocoann ^x coann (r> i n\ 

IT - n ' 1 j 

versus the mass ratio of the lightest coannihilating particle and the neutralino in our sample 
of models: the coding in relic density shows that, for cosmologically interesting cases, only 
mass differences up to about 40% are important, hence the cut at the mass ratio 1.5 (vertical 
dashed line in the figure) is sufficiently conservative. If a 1% accuracy is required even for 
models that are cosmologically disfavoured, one would need to raise the cut to about 
1.7: the coannihilation effects one picks in this way, however, correspond to the Xi ~ X? 
annihilation cross section being very suppressed and, even adding coannihilation terms, 
the relic abundance still remains r2 x /i 2 S> 1. In Fig. |4| we show the effective annihilation 
cross section for such an example (model C in Table [2| in Appendix |C]). In this case stop 
coannihilations are important even though the stop mass is about 50% higher than the 
neutralino mass, but still the relic density, which is £l x h 2 = 73.2 when coannihilations are 
neglected, is shifted down to just f^/i 2 = 19.7. The criterion with selection on mass we 



- 10 - 



a 



1000000% 



Cj 100000% 
10000% 
1000% 
100% 
10% 

1% 

1%o 



-2 

XXXXXXXXX> 

xxxxxxxxxx 
xxxxxxxxxx 





J. Edsjd, M. Schelke, P. Ullio and P. Gondolo, 2003 



AO, positive 



-_ 

xxxx 
xxxx 



§ XX 
**XX 
**XX 
+ + + 
+ 4-4-4- 
+ + + + + 
+ + + + + 
+ + + + + 
+ + + + 4 



XXXXX> 
XXXXX) 

xxxxx> 



XXXXX> 
XXXXX) 
XX**X> 
+ + H 
+ + H 



+ + + + 

xx**4 
xxx*4 

XX* + * 
+ + + 4 
+ + + 

x***4 

X** + + 
4- + 4 
+ + 4 



++++++++++++++++++ 



+++++++++ 

- ++++++++ 
++++*** 



H»€»0 8! X 

B©$OOOX 

h + + O X 

h + 

b©ooooooooox 
e«oo»o o oooc 
&©ooooo oo o o© 

a +_o__ »g__' 

r^F JU GO O OO 

h©000+0 OO O O 

h©eeo o oo oo 

h©*ooo ooo ooooooooo 





<^ir) ( 
0.1 < (Oh 2 ) 



<0.1 



. <0.3 



4- 4- 4- 4-4- 4- 4- 4- 4-4 



oo o o ooo o 

DOOOO OOOO O 
ZJOOO OO OOOO 

bo©oooo oo oo 
Boeooooo o oo 



S'oioeiis 

0000030$ 



4-4-4- 4-4- ******* 



oooooc 
ooookc 
oooooc 

I ! ! ' 



BBOOOOO OOOOl 
«»SOOOOO o < 
***4- jOOOO Ol 
4- + 4- 4- -ft * + 



+ ++ + 4- + + + + 4- + + + + 4- + + + 



4-+ 4- 4- 4- 4- + 4- 4- 4- + 4- 4- + 4- 4- 4- 4- 4-4-4[ 



"1 | + + + * + * r 

z 000000*++*++ +*+«000®0000000 OOOOOS 
4- + + + + + + 4-4- OOOOOO OCj 

- ooooooa©e©®4- oooooooooooo 

mo/ ± 4-4- 4-4-** 

- I /O = + + + * + 

= oo»o©®©+ 



AO. negative 



4-4-4-*4-X 



-10% 



-100% 



1.1 



1.2 1.3 



1.4 



1.5 



1.6 



1.7 



^lightest coannihilation particle ^ 



Figure 3: The relative change in relic density due to the inclusion of coannihilations, Afi/fi = 
(f^no coann — ^coann)/^coann, versus the mass ratio between the lightest coannihilating particle and 
the lightest neutralino. To gain computational speed we avoid including unnecessary coannihilation 
processes in the calculation by imposing a cut at 1.5m x . The models included in the figure are 
from some general mSUGRA scans with m <E [0,2500] GeV, m 1/2 G [0,2500] GeV, tan/3 S [5,50], 
A e [-5000, 5000] GeV and sign(^) = ±. All models in Figs. [l0|-|l4] are also included here. 



give is the easiest to implement in the numerical calculation, but clearly depends on the 
strength of the interactions we are dealing with; we can say that it is perfectly safe in the 
scheme we are dealing with, but need not be valid in other schemes. 

A final check we perform is the following: sometimes in the literature only the next-to- 
lightest (NLSP) sparticle has been included in the coannihilation calculation and one may 
question whether this can be considered a fair approximation. There is reason to believe 
that also heavier sparticles might change the relic LSP density by a non-negligible amount 
if they are either close in mass to the NLSP or if they have larger coupling strength than 
the NLSP (as e.g. nas compared to Whether or not it is a good approximation 
to include only NLSP coannihilations should therefore be checked case by case. Even 
restricting to the cases of coannihilating particles with comparable couplings, such as staus, 
smuons and selectrons, it is not straightforward to provide a firm criterion telling when to 
include all of them or just the NLSP (i.e. the f in this case) in the calculation. In practice, 
we find regions in the mSUGRA parameter space when smuon and selectron contributions 
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Figure 4: The effective annihilation cross section as a function of p e g for an example model (model 
C in Table ^ in Appendix ^) where stop coannihilations are important even though the lightest 
stop is about 50% heavier than the lightest neutralino. 
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Figure 5: The effective annihilation cross section versus p c g for an example (model D in Table g 
in Appendix ^|) where inclusion of all three lightest sleptons is important to get the correct relic 
density. 



give corrections to the relic density of the order of 10% or higher. An example is given 
in Fig. H (model D in Table || in Appendix P). In this example, the lightest smuon and 
selectron are almost degenerate in mass with both the lightest stau and lightest neutralino. 
Including in the calculation the coannihilation with these three sfermions result in a relic 
density Q x h = 0.109. If instead we only include ypf and ff coannihilations, the relic 
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density would be 0.128, i.e. we would be ~ 18% off from the correct value. Hence, as this 
example has shown, it can be important to include the coannihilations of more than one 
particle. To be on the safe side, we always include all particles with masses up to 1.5m x . 



4. The role of coannihilations in the mSUGRA framework 

The mSUGRA framework has been extensively discussed in the literature. We sum- 
marize here the main features that are relevant for understanding the results of the relic 
abundance computation. In line with most previous analyses, we sample the 5-dimensional 
mSUGRA parameter space choosing a few values of tan (5 and Aq, and slices along the 
m 1 / 2 ,mo planes for both sign(n). 

Consider first the case Aq = 0. Two regimes with cosmologically interesting relic 
abundances have been identified: The region mo < Tn\/2 and the region mo 3> fn\/2- In 
the first region, the lightest neutralino is a quasi pure bino with mass set essentially by m^ 2 
alone; the parameter mo sets the sfermion mass scale, with the slepton sector lighter than 
the squark sector and with the lightest stau always being the lightest sfermion, possibly 
lighter than the lightest neutralino if mo <C m\i%. In the second regime, the relevant region 
is a narrow band, sometimes dubbed the "focus point" region |43[| , close to the region where 
there is no radiative electro-weak symmetry breaking: in such a band the parameter n is 
driven to small values and forces a mixing between the gaugino and Higgsino sectors. As 
a consequence the lightest neutralino may contain a large Higgsino fraction and the next- 
to-lightest supersymmetric particle is a chargino. Large values of Aq can introduce a third 
regime in the intermediate mo range: off-diagonal entries in the stop mass matrix can 
become sufficiently large to drive the lightest stop to masses smaller than the mass of the 
lightest stau or even the mass of the lightest neutralino. On top of these generic trends, 
details are sensitive to the value of tan/3 and the sign(fi). 

In each of the regimes above, the mass splitting between the LSP, which we require 
to be the lightest neutralino, and the next-to-lightest SUSY particle can be small enough 
for coannihilation effects to become important. We label the three regions as the slepton 
coannihilation region, the chargino coannihilation region and the stop coannihilation region, 
and describe each of them separately. 

4.1 Slepton coannihilations 

We consider first the case Aq = and mo < n^-i/2- The slepton coannihilation region 
was recognised in Ref. [|l4| and has been the focus of several recent studies, including 
|15, 13, 21, 44]. As an example, we take tan/3 = 10 and study the mi/2 — mo plane for 



both signs of /i. In the top panels of Fig. y, we plot isolevel curves for the neutralino relic 
abundance, including coannihilation effects, for 8 different values of Q x h 2 starting from 
£l x h 2 = 0.3 and decreasing down to £l x h? = 0.025. The latter value corresponds to the 
case in which neutralinos would be a subdominant dark matter component in the Universe 
but could still account for a major part of the dark matter in galaxies. The shaded area on 
the left in each panel is excluded by accelerator constraints, while the shaded area towards 
the bottom right corner in each panel is removed because in this region the LSP is the 
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Figure 6: Results for tan/3 = 10 and Ao — 0. The isolevel curves for the relic density Q, x h 2 are 
shown in the top panels. In the bottom panels, curves indicate how big the error on the relic density 
would be if coannihilations were not included. The mass splitting between the lightest neutralino 
and the lightest stau is also indicated. 

lightest stau rather than the lightest neutralino: its upper bound marks the line along 
which the (bino-like) neutralino and the lightest stau have equal mass. 

We can give a schematic interpretation of the results displayed starting with the isolevel 
curves on the top left corner of each panel, where all isolevel curves converge to a narrow 
band. There, the model has a relatively heavy sfermion sector, and the lightest neutralino 
mass is just a few GeV larger than half the Z° boson mass. The bino pair annihilation rate 
into fermions is dominated by the diagram with Z° in the s-channel at energies just slightly 
displaced from the Z° resonance: this resonant annihilation leads to acceptable values of 
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Q x h 2 in a narrow band. Thus Q x h 2 is very sensitive to the parameter mi/2 as the value of 
the bino mass has to be fine-tuned in such a way that the thermally averaged cross section 
picks the right portion of the resonance; for larger m x , (av) drops rapidly and Q x h 2 becomes 
large, closer to the resonance (av) becomes exceedingly large and VL x h 2 very small. When 
we follow the isolevel curves from the top left down to smaller values of mo, sfermion 
masses decrease and the amplitude for neutralino annihilations into fermions mediated 
by sfermion exchange in t- and u-channels increase, eventually becoming dominant. The 
largest contribution to (av) is given by the t~t + final state which is mediated by the lightest 
sfermion, in this case the lightest stau. Such an increase in the cross section when moving 
to lower mo can be compensated by increasing the bino mass, i.e. by shifting to larger 
This explains the trend of the isolevel curves starting to align along the diagonal in the 
figure. Further down this diagonal the mass of the lightest stau is more or less constant, 
but the mass of the heaviest stau increases. This increased mass splitting between the 
staus causes the cross section to increase, but this increase is again compensated by a 
larger neutralino mass. Before reaching the lower right corner where the lightest stau is 
the LSP, coannihilation effects take over: (a e gv) becomes rapidly dominated by neutralino- 
slepton and slepton-slepton contributions. The coannihilation cross section decreases with 
increasing initial state masses, an effect which can be compensated by decreasing the mass 
splitting between the LSP and the stau. (This effect is further discussed in connection 
with Fig. ||[) Hence, the isolevel curves bend almost parallel to the bound of the excluded 
region, confined to a band which gets progressively narrower towards larger values of mo 
and m-L/2- As can be seen by comparing the top left with the top right panel, the flip in 
sign(fi) does not alter this overall picture but just slightly displaces the position of the 
isolevel curves in the mi/ 2 — tuq plane. 

The region where coannihilation effects become important is highlighted in the bottom 
panels of Fig. ^. There dashed lines show the isolevel curves for the relative difference in relic 
abundance computed neglecting and including coannihilations, Afl/fi = (f2 X) no coann — 
^x, coann)/^x> coann- The values we display span from 1% to 500%. Isolevel curves for the 
relative mass splitting between the lightest stau and the lightest neutralino are shown as 
well and, as expected, the correlation with the isolevel curves of Afi/fl is evident. 

In Fig. [?] we repeat the exercise for tan (5 = 30 and the picture we see is analogous to 
what we found in the tan/3 = 10 case. The shift of the isolevel curves to slightly larger 
niQ values is mainly due to the fact that, comparing corresponding points in the plane, 
sleptons are driven to slightly lighter masses at larger tan (3. 

From Figs. |6| and |7[ as well as from other sample checks performed for other values 
of tan/3, we can infer that, as a rule of thumb, in the Aq = case a 25% mass splitting 
between the bino and the lightest stau is approximately the borderline below which slepton 
coannihilations have to be included for an estimate of £l x with a 1% accuracy level. 

We have briefly mentioned the reason why the stau coannihilation region takes the form 
of a tail extending to large values of mj^. When m 1( / 2 , and therefore m x , is increased, it is 
necessary to increase the effect of coannihilations, and therefore to lower the relative mass 
splitting between the LSP and the NLSP, if the relic density should not increase. Fig. || gives 
a detailed illustration of this effect. We select here two models with Aq = 0, tan/3 = 10, 



- 15 - 





Figure 7: Results for tan/3 = 30 and Aq = 0. The isolevel curves for the relic density Q, x h 2 are 
shown in the top panels. In the bottom panels, curves indicate how big the error on the relic density 
would be if coannihilations were not included. The mass splitting between the lightest neutralino 
and the lightest stau is also indicated. 



sign(fj) positive, and the m\/2i mo pair chosen in such a way that the relic abundance 
for both models is Q x h? = 0.115, but AQ/Q is 100% for one model (left panel, model 
E in Table and 1000% for the other (right panel, model F in Table |2|). Analogously 
to some of the figures in Section |3.2| , we plot W e g/4:E^ S , as well as the weight function 
k computed at the freeze out temperature. Going from the left to the right panel, m x 
increases from 138.5 GeV to 371.1 GeV, the lightest stau from 148.0 GeV to 371.8 GeV. 
The mass splitting consequently goes from 6.8% to 0.21%. The increase in neutralino mass 
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Figure 8: The effective annihilation cross section W c g / ^E^ s versus p c s for two sample cases (model 
E and F in Tabic | in Appendix |C]) with tan/3 = 10 and Q, x h 2 = 0.115. The fi eure illustrates how 
coannihilation effects become more important the further out in the coannhilation tail we get. See 
text for a discussion. 

causes a sharp drop in the Xi~Xi annihilation cross section 7 ; on the other hand such a 
drop can be compensated in the thermally averaged cross section by increasing the role 
of coannihilating particles. Even though their annihilation cross sections have also been 
slightly reduced by the mass shift, we have forced their contribution to the thermal average 
to be larger by shifting them deeply in the region where the weight function k is large. 
More precisely, we have required the coannihilation threshold to move to lower p c g, i.e. to 
have a lower mass splitting. 8 With these changes in the mass splitting, we have found two 
models with about the same thermally averaged cross sections, and then, with the rule of 
thumb Q x h 2 ~ 10 -27 cm 3 s -1 /(<7 e ffu) one can expect similar relic densities. In fact, for the 
two cases discussed here we have chosen the parameters to give exactly the same neutralino 
relic density £l x b? = 0.115 from the full calculation. 

7 As W c ff/4_E^ ff at p c ff = is the neutralino annihilation rate at zero temperature, which is the relevant 
quantity for indirect dark matter detection, we find that, roughly speaking, it is a factor of 140 harder to 
detect indirectly the model on the right, a factor of 20 due to the reduced value of the cross section, plus a 
factor of 7 due to the square of the ratio of m x 's (which takes into account the neutralino number density 
scaling for a fixed dark matter mass density). 

8 Note that if we had kept the mass splitting fixed and just shifted the mass scale, we would have found 
an equal shift in the position of the coannihilation thresholds on the p e ff scale. On the other hand, the 
"width" in p e ff of the weight function k would have increased by about the same factor since the freeze-out 
temperature is higher. Hence, the net result would have been that the right-hand panel would have looked 
very similar to the left-hand panel, except for roughly an overall stretch along the p c g scale, and for all cross 
sections being lower. We remind the reader that we have fixed the range of the scale on the top of each 
panel (the value on the scale marks the location of the momentum p^ defined as «(p^)/K(p"Jf x ) = 10", 
where p"ff x is the momentum corresponding to the maximum of the weight function k) and derived the 
corresponding range in p e s. 
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Figure 9: For a set of values of £l x h 2 , we show the importance, Af2/0, of including slepton 
coannihilations, as well as the mass splitting between lightest neutralino and lightest stau, both as 
functions of the neutralino mass. We have here chosen tan/3 = 30, fi > 0, A n = and low m (the 
latter to be in the stau coannihilation region). In the right-hand panel, we also indicate by arrows 
where the upper limit on the mass would be if coannihilations were not included (for these limits 
the vertical axis is meaningless). The shift to higher masses when coannihilations are included is 
clearly seen. 



From the above discussion it should also be evident that we cannot play the same game 
to arbitrarily large neutralino masses. At some point one reaches the edge at which the role 
of coannihilating particles cannot be further strengthened and one finds an upper bound on 
the neutralino mass (for that specific choice of tan/3, signal) and Q x h 2 ). The introduction 
of slepton coannihilations have extended the cosmologically allowed mass interval of the 
bino-like LSP in the mSUGRA context. This fact has already been pointed out by several 



authors. In the work of some authors [20, [44[] it seems that the neutralino mass might be 



unbounded from above, but the majority [21, 14, 15] found that there was still an upper 
limit to the allowed mass. Given the high precision in the relic abundance calculation we 
present in this analysis, we are able to make a qualitative analysis of the high m x region, 
and indeed we find upper limits to the neutralino mass. We will show our results in the 
case tan/3 = 30 and positive \i considered in Fig. |?j. Keeping in mind that we are in 
the mo < mi/2 regime, there is a one to one correspondence between the pair (m 1( / 2 ,mo) 
and the pair (m x , (jrif — m x )/m x ) or (m x ,Ail/fl). We therefore show first in Fig. ||a the 
difference in relic density, (r2 XjnoC oann — ^x>coann)/^x>coanru versus the neutralino mass for 
a few values of Q x h 2 . We clearly see the importance of stau coannihilations in the high 
mass region. In Fig. we instead show the relative neutralino-stau mass splitting (i.e. 
Am = {rrif — m x )/m x ) as a function of m x . Also shown in the figure are the upper limits 
on m x for the case where Q x h 2 is computed ignoring coannihilation effects; the shift to 



- 18 - 




Figure 10: The relic density contours (solid lines) for models in the focus point region; tan/3 = 30 
and Aq — 0. In a) /i > and in b) /i < 0. The kinematic chargino mass limit of 104 GeV and the 
W + W~ and tt thresholds are indicated. 

much larger values, when coannihilation effects are included, is evident, as well as the fact 
that we do find a new maximum value of m x . The results for negative \x are very similar, 
while those for tan/3 = 10 are analogous, but show slightly more stringent upper bounds 
on the neutralino mass. 

4.2 Chargino coannihilations 

In the focus point region, the value of the soft mass parameter (at the electro-weak scale) 
for the Higgs doublet that couples to up-type quarks, ulh u , is naturally of the electro- weak 
scale, regardless of mo [|43|| . As a consequence, the parameter [i is forced to be light, and 
can be at the level of the gaugino mass parameter mi/2 or even lower. This implies that the 
neutralino LSP may have a large Higgsino fraction and be nearly degenerate in mass with 
the lightest chargino and the next-to-lightest neutralino. Especially at higher mo-values, 
the Higgsino fraction can be very large, close to one. Hence, in this high mo focus point 
region, chargino (and neutralino) coannihilations are expected to be important. Chargino 



coannihilations have been extensively studied in the generic MSSM context [45, [5], 10], 
but have been rarely stressed in the mSUGRA framework (although they are included in 
some recent analyses, e.g. [20, |47|]). 



In Fig. [lO| we show the lower part of the focus point region for tan/3 = 30 and A§ = 0. 
The top-left corners of these figures are excluded due to no radiative electro- weak symmetry 
breaking, but close to that region, we see the focus point region emerge. In this region, 
the Higgsino fraction is usually small, but non-negligible, and the same is true for the 
effect of chargino coannihilations. This is the part of the focus point region most often 
discussed in the literature. However, if we continue to higher masses, we get a band of 
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Figure 11: The relic density contours (solid lines) for high mass models where chargino coanni- 
hilations are important; tan/3 = 30, Aq — and /i > 0. In a) coannihilations are not included, 
whereas they are included in b). Neutralino mass contours are shown with dashed lines. 




Figure 12: We here show the effect of coannihilations for models where chargino coannihilations 
are important; tan/3 — 30, /i > and Aq — 0. We show in a) Af2/f2 = (fl no coallI1 — f2 C oann)/^coann 
versus the neutralino mass and in b) the mass splitting between the lightest chargino and the 
lightest neutralino versus the neutralino mass. 



cosmologically interesting relic densities where the Higgsino fraction increases as we go 
up in mass (at the highest masses, it is close to 1). In this case, coannihilations with 
the lightest chargino (and the next-to-lightest neutralino (s)) occur and are important. In 
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Fig. [Tl] a we show the relic density isolevel curves without including coannihilations and in 
b) with coannihilations included in this high-mass region. (Note the scale mo-2m^ 2 on the 
y-axis, which is chosen to clearly show the coannihilation region.) We clearly see the effect 
of coannihilations pushing the cosmologically allowed region to higher masses. Note that in 
the focus point region of the parameter space, there is no longer a simple relation between 
the neutralino mass and mi/2- Some neutralino mass contours are therefore indicated in 
the figures. 



As shown in Fig. 12, we find that (just as for the MSSM case [|1C| 1) chargino coanni- 
hilations are non-negligible below the W mass where the dominant neutralino-neutralino 
annihilation channel is annihilation into fermion anti-fermion pairs via s-channel Z-boson 
exchange. The coupling for a Higgsino-like neutralino to the Z-boson is very suppressed 
though, whereas that of charginos is not and coannihilations could then give a big effect. 
However, in the mSUGRA framework, the neutralinos below the W mass are rather mixed 
than Higgsino-like and the effect is not as dramatic as in the MSSM. Above the W mass, 
where annihilation into W + W~ dominates for neutralino-neutralino annihilation (annihila- 
tion into Z°Z° is also significant above the Z mass), the effect of coannihilations is smaller 
since the annihilation rate into W + W~ is comparable for neutralinos and charginos. Still, 
annihilation into fermion anti-fermion pairs is not suppressed for chargino-chargino annihi- 
lation and thus the chargino-chargino annihilation cross section is typically slightly higher 
than the neutralino-neutralino one. Hence, coannihilations do change the relic density even 
for higher masses, and the further up in mass (and down in mass splitting) we go, the more 
important they are. In Fig. U a we plot AO/O versus the neutralino mass. We see that 
AO/17 is of the order of 1% for models with m x just below the W mass and drops when 
we get above m\y- Then, AO/O slowly increases as we go up in mass with a small drop at 
m x = 175 GeV, where the neutralino annihilation into ti becomes significant, thereby re- 



ducing the importance of coannihilations somewhat. In Fig. 12b we show the relative mass 
splitting between the lightest chargino and the lightest neutralino versus the neutralino 
mass. We here clearly see the same effect as we have seen with slepton coannihilations, 
i.e. as we go up in mass, we need to make the mass splitting smaller to maintain the same 
value of the relic density. As for slepton coannihilations, we cannot continue this game to 



arbitrarily high masses and as we see in Fig. |12p, we get upper limits on the neutralino 
mass for a given relic density. The corresponding curves for /j, < are very similar and we 
do not show them separately. 

4.3 Stop coannihilations 

The new version of DarkSUSY offers the possibility to include all squark coannihilations 
in the relic density calculations, but only stop coannihilations have proven to affect the 
result for the mSUGRA framework. In this context, the lightest stop is always the lightest 
squark; its mass is usually much larger than the LSP mass, unless off-diagonal entries in 
the stop mass matrix become sufficiently large to drive the lightest stop mass to small 
values: this can happen when \Aq\ is large. If the mass of the lightest stop is close to 
the neutralino mass, stop coannihilation effects become important in the neutralino relic 
density calculation [12|, 16, 45]. In Fig. 13 we show the relic density in the m^j — mo-plane 
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Figure 13: Isolevel curves for the relic density a) without coannihilations and b) with coannihi- 
lations for an example of mSUGRA parameters where stop coannihilations are important. In the 
hatched region to the left, t\ is the LSP (except for a narrow band at the lower part of the right 
edge where the neutralino is the LSP, but where the H® mass constraint is not fulfilled). In the 
hatched region at the bottom right, fi is the LSP. In the almost vertical band of interesting relic 
densities, a light t\ is important in two respects, both by boosting the Xi ~ Xi annihilation above 
the tt threshold (as seen in a)) and by coannihilations (as seen in b)). The most prominent effect 
of the ti coannihilations is the narrow band for neutralino masses less than m(t) ~ 174 GeV. 



for tan /3 = 10, Aq = —2750 GeV and /i > 0. The hatched region is excluded mainly due to 
the neutralino not being the LSP here. The lower right part has the fi as the LSP, while ti 
is the LSP in most of the 'bump' to the left. A narrow band close to the lower part of the 
right edge of this 'bump', i.e. to the left of the Qh 2 curves, is excluded because of the mass 
limit on H®. In other words, the H® mass limit cuts away some part of the parameter space 
with the smallest mass splitting between t\ and the LSP. In Fig. 13a we show the relic 
density isolevel curves without including coannihilation processes. We have also included 
the isolevel curves for the neutralino mass. It is clearly seen that the neutralino-neutralino 
annihilation cross section dramatically increases at the ti threshold. The relic density then 
decreases to cosmologically interesting values in a band were the mass splitting between 
the t\ and the neutralino is small. It is the t-channel exchange of the light stop that 
boosts the neutralino annihilation into ti. In Fig. 13b we show the results of including 
coannihilation processes. Three effects are seen. First of all, the stop coannihilations open 
a narrow band of cosmologically interesting density for neutralinos below the top mass. 
Secondly, in the almost vertical band above m(xi) = m (t), that was already present in 
Fig. |l3|a, the relic density is decreased further due to t± coannihilations. Finally, we also 
find interesting relic densities for high m^ 2 and low mo where fi-coannihilations lower 
the relic density. The same three features were also found by Ellis, Olive and Santoso, 
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Figure 14: For the same set of models as in Fig. |l3], we plot ACl/Cl = (O no coanu — f2 C oann)/^coann 
versus the neutralino mass. We clearly see the effect that below mt, we only get interesting relic 
densities if coannihilations with t are included. When we get close to mt, though, even XiXi~ 
annihilation gets efficient and the relative effect of including coannihilations go down. When coan- 
nihilations with f get important at higher masses, the relative difference in Cl x goes up again. The 
small dip at m x ~ 280 GeV corresponds to the low-mo corner where the curves in Fig. [l3| bend 
sharply and arises because XiXi-annihilation to t + t~ gets significant here. 



1 16] Fig. 6d. Beside the trivial sign difference in the convention used for Aq, we also have 
chosen the numerical value of Aq slightly lower than done in [16], in order to obtain a 
figure as similar to theirs as possible. The possible explanation for the required change of 
\Aq\ is that we use different RGE codes, and a given |^4q| therefore results in different low 
energy masses. This is not only true for the neutralino and sfermion masses but also for 
the Higgs masses. As mentioned, part of the edge of the excluded region is set by the H% 
mass limit. The excluded region would be even larger if we applied the Standard Model 
Higgs limit (~ 114 GeV) as done in [16] instead of the value 771(1/2) ~ 92 GeV that applies 
for tan (3 = 10 |32 



To see more clearly the effect of coannihilations, in Fig. |14J we show the difference in 



relic density, (CI 



X, no coann 



Cl 



X, coann 



X, coann; 



versus the neutralino mass along the isolevel 
curves shown in Fig. 13. We see that in the coannihilation region in the upper left corner 
of Fig. 13b, ti-coannihilations are indeed very important changing the relic density by a 
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factor of several 100. When we get to the ti threshold, the effect of coannihilations decreases 
drastically since the X1X1 annihilation cross section is high. When the neutralino mass is 
increased further, the importance of coannihilations again increases and above ~ 275 GeV 
the f-coannihilations start getting important. 

5. Conclusions 

We have presented a novel tool for calculating the neutralino relic abundances with an 
estimated precision of 1% or better, assuming masses, widths and couplings of particles in 
the MSSM are given. We allow for the most generic coannihilation effect in the framework 
of the MSSM, applying at the same time the state of the art technique to trace the freeze- 
out of a species in the early Universe. The code for numerical computations will be publicly 
available in the near future together with a new extended version of the DarkSUSY package 

As a first example, in this paper we have discussed results valid in the mSUGRA setup, 
the rather restrictive framework most recent analyses of coannihilations have focussed on, 
but, at the same time, general enough to illustrate most of the effects that can arise 
in the MSSM. For a selection of mSUGRA models we have given the reader a visual 
guide into general trends and their exceptions, such as the tendency of coannihilation 
effects to lower the neutralino relic abundances and the cases in which the opposite can 
happen. We have then performed a broad study of the neutralino relic density in the 
mSUGRA parameter space. The features we have spotted are in agreement with the picture 
emerging from previous analyses; in particular we have discussed the cases of neutralino 
coannihilations with sleptons (most notably stau leptons), with stop squarks, and with 
charginos. Especially the chargino coannihilations have been treated in more detail than 
before. 

The accuracy of the calculation we performed allowed us to go into great details in 
the cases presented. Novel features we have discussed include, e.g.: (i) a rule of thumb 
to discriminate the case when slepton coannihilations are relevant: 25% mass splitting 
between lightest stau and lightest neutralino for a 1% accuracy on the cosmologically 
allowed neutralino relic density; (ii) a shift to heavier masses of the cosmological bound on 
the neutralino mass in both the slepton and chargino coannihilation cases: in the tan f3 = 30 
case we studied, we found 565 GeV as an upper limit on the neutralino mass from the 
cosmological bound ^ x /i 2 < 0.2 in the regime mo < mi/2, where slepton coannihilations 
are important, and 1500 GeV in the regime mo 3> mi^, where chargino coannihilations 
occur; (iii) a correlation between mass splittings, differences in VL x h? and neutralino masses 
in all of the cases presented. 

In an upcoming paper, we will investigate the effects of the coannihilation channels in 
a more general MSSM context. 
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A. Included coannihilations 

Here we list all 2 — > 2 tree-level co-annihilation processes with sfermions, charginos and 
neutralinos. All the processes are included in the DarkSUSY code, but not all of them have 
been used in the calculation we report on. As mentioned earlier, we have found that an 
accuracy of 1% on the relic neutralino density in the cosmologically interesting interval is 
obtained by including coannihilations for particles lighter than 1.5m x . Consequently, some 
sparticles never get included in the initial state. For the sparticles that satisfy the mass 
difference criterium, we have included all coannihilation processes and for each of these, all 
the exchange channels. Furthermore, the DarkSUSY relic density code always includes the 
one loop neutralino annihilation into 77, 7 Z and gg. 

It should be noted that we have not included all flavour-changing charged current dia- 
grams. The DarkSUSY vertex code for the charged current couplings is written in a general 
form that includes all possible flavour-changing (and flavour-conserving) vertices. The 
flavour-conserving couplings are much larger than the flavour-changing. For the sfermion 
coannihilations with charged currents we only take the flavour-conserving contributions, 
while for the chargino coannihilations we include the flavour-changing contributions as well. 
In a future version of DarkSUSY, we may as well include the flavour-changing processes for 
the sfermion coannihilations, even if they are not expected to be important. 

We have used the notation / for sfermions and / for fermions. Whenever the isospin 
of the sfermion/fermion is important, it is indicated by an index u (T3 = 1/2) or d 
(T3 = —1/2). The sfermions have an additional mass eigenstate index, that can take 
the values 1 and 2 (except for the sneutrinos which only have one mass eigenstate). A 
further complication to the notation is when the sfermions and fermions in initial, final 
and exchange state can belong to different families. Primes will be used to indicate when 
we have this freedom to choose the flavour. So, e.g. f u and f u will belong to the same 
family while f u and f' u can belong to the same or to different families. Note that the colour 
index of (s)quarks as well as gluons (g) and gluinos (g) is suppressed. 

Besides the sfermions we also have neutralinos and charginos in the initial states. The 
notation used for these are the following. The neutralinos are denoted by x? with the index 
running from 1 to 4. The charginos are similarly denoted Xj with the index taking the 
values 1 and 2. 

In the table in this appendix, a common notation is introduced for gauge and Higgs 
bosons in the final state. We denote these with B with an upper index indicating the 
electric charge. So B° means H® , H® , #3 , Z, 7 and g while is and . We will 
use additional lower indices m and n when we have more than one boson in the final state. 
Thus indicating that the bosons can be either different or identical. Note that the case of 
two different bosons also includes final states with one gauge boson and one higgs boson. 

The table has been made very general. This means that when a set of initial and 
final state (s)particles have been specified, the given process might not run through all the 
exchange channels listed for the generic process. Exceptions occur whenever an exchange 
(s)particle does not couple to the specific choice of initial and/or final state. As an example 
we see that since the photon does not couple to neutral (s)particles, none of the exchange 
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channels listed for the generic process fi + x^ B° + f actually exist for the specific 
process v + x° — ► 7 + v. All these exceptions can be found in the extended tables in Ref. 



1 49]. Also note that the list of processes is not complete with respect to trivial charge 
conjugation. For each process of nonvanishing total electric charge in the initial state there 
exist another process which is obtained by charge conjugation. 
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Table 1: Included co-annihilation processes through s— , t— , u— channels and four-point interactions 
(p). For the /t^/J processes the corresponding process for up-type sfermions can be obtained by 



interchanging the u and d indices. 
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B. A note about internal degrees of freedom 



Here we describe a technical detail in the calculation, which we find useful to specify. If 
we look at Eqs. ( |3.5| ) and ([O]) we see that we have a freedom on how to treat particles 
degenerate in mass. For example, the charginos, x^i (where the mass eigenstate index is 
implicit) can be treated either as (a) two separate species x + an d X~, each with internal 
degrees of freedom g x + = g x ~ = 2, or (b) a single species x^ with g x ± = 4 internal degrees 
of freedom. Of course the two views are equivalent, we just have to be careful to include 
the giS consistently whichever view we take. In a), we have the advantage that all the W« 
that enter into Eq. ( |3.5| ) enter as they are, i.e. without any correction factors for the degrees 
of freedom. On the other hand we get many terms in the sum that are identical and we 
need some book-keeping machinery to avoid calculating identical terms more than once. 
On the other hand, with option b), the sum over Wij in Eq. ( |3.5| ) is much simpler only 
containing terms that are not identical (except for the trivial identity W%j = Wji which is 
easily taken care of). However, the individual Wij will be some linear combinations of the 
more basic Wij entering in option a), where the coefficients have to be calculated for each 
specific type of initial condition. 

We have chosen to work with option b) since this most easily gives an efficient numerical 
code. Denoting the W^s in option b) with a prime, we can derive 4S] the following relations 
between the two different views, 
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Where / denotes a generic sfermion. We have neglected in t 



averaging for squarks, i.e. W~, ~, = \\ Y^ a b=i 



(B.l) 

lis listing the trivial colour 
and similarly for the other 



cases. 
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C. Example models 

Here we list the parameters and some properties of the example models considered in some 
of the figures. 
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Table 2: Model parameters and some properties of the example models discussed in the text and 
figures. The sparticle masses are calculated with ISASUGRA 7.64. 
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